Introduction to Spatial Analysis in R

Material adapted from Grant R. McDermott’s notes

NOTE: the examples in this notebook focus only on vector-based spatial analysis. For information on raster-based spatial analysis, check some of the resources from Data Carpentry:

https://datacarpentry.org/organization-geospatial/

https://datacarpentry.org/r-raster-vector-geospatial/

Some Tools for Spatial Analysis in R

External libraries

(requirements vary by OS)

R provides bindings to powerful open-source GIS libraries. These include the - Geospatial Data Abstraction Library (GDAL) - Interface to Geometry Engine Open Source (GEOS) API suite - Access to projection and transformation operations from the PROJ.4 library.

R packages

• Some new R packages for this introductory set of examples: sf, lwgeom, maps, mapdata, spData, tigris, tidycensus, tmap, tmaptools

Note that you need ggplot2 version 3.0.0 or above for the required plotting support of sf objects.

library(tidyverse)
library(tidycensus)
library(sf)
library(tmap)
library(maps)
library(mapdata)
library(tigris)

Census API key

We will access some data from the US Census Bureau through the tidycensus package. This will require a Census API key, which you can request here. Once that is done, you can set it using the tidycensus::census_api_key() function. You can use the "install = T" option to save your key.

Introduction: Coordinate Reference System (CRS) and map projections

Additional References: - Overview of Coordinate Reference Systems (CRS) in R

• Spatial Data Science with R

• Geocomputation with R

Spatial data (like all coordinate-based systems) only make sense relative to some fixed point. That fixed point is what the Coordinate Reference Systems, or CRS, is trying to set. In R, we define the CRS with a “proj4string”, which is based on the syntax of the PROJ.4 geospatial library.

In R, the notation used to describe the CRS is proj4string from the PROJ.4 library. It looks like this:

+init=epsg:4121 +proj=longlat +ellps=GRS80
+datum=GGRS87 +no_defs +towgs84=-199.87,74.79,246.62

Whenever we try to plot (some part of) the earth on a map, we are effectively trying to project a 3-D object onto a 2-D surface. This will necessarily create some kind of distortion. Different types of map projections limit distortions for some parts of the world at the expense of others. Distortions are particularly acute in the extreme latitudes for some projections, so try to choose according to the specific needs of your research question. For example, take a look at how badly the (in)famous Mercator projection distorts parts of the global map (source)

Animating the Mercator projection to the true size of each country in relation to all the others.

Focusing on a single country helps to see effect best.#dataviz #maps #GIS #projectionmapping #mapping pic.twitter.com/clpCiluS1z

— Neil Kaye (@neilrkaye) October 12, 2018

You can also “orient” your projection to a specific latitude and longitude using the PROJ.4 syntax.

Simple Features and the sf package

R has long provided excellent support for spatial analysis and plotting (primarily through the sp, rgdal, rgeos, and raster packages). However, until recently, the complex structure of spatial data needed a set of equally complex spatial objects in R. The sf package provides simple features access for R. The “sf” stands for simple features, which is a simple standard for representing the spatial geometries of real-world objects on a computer. See the documentation and excellent sf vignettes for more details.

These objects — i.e. “features” — could include a tree, a building, a country’s border, or the entire globe. The point is that they are characterized by a common set of rules, defining everything from how they are stored on our computer to which geometrical operations can be applied them. Of greater importance to us here, however, is the fact that sf represents these features in R as data frames. This means that all of our favorite tidyverse operations can be applied to spatial data.

Reading in spatial data

Most of the functions in the sf package start with the prefix st_. This stands for spatial and temporal, and a basic command of this package is easy enough once you remember that you are probably looking for st_SOMETHING().

As a quick demonstration we will read in the Florida counties shapefile. We use the st_read() command and the sf package will handle all the heavy lifting behind the scenes.

The shapefile was downloaded from: http://geodata.myflorida.com/datasets/swfwmd::florida-counties

# library(sf) ## Already loaded
## Location of our shapefile (here: bundled together with the sf package)
loc_counties <- "./Florida_Counties/Florida_Counties.shp"
## Read the shapefile into R
fla_counties <- st_read(loc_counties, quiet = TRUE)

Simple Features as data frames

Print out the fla object that we just created and take a look at its structure.

fla_counties

Simple feature collection with 67 features and 7 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -87.62601 ymin: 24.54522 xmax: -80.03095 ymax: 30.99702
CRS:            4326
First 10 features:
   OBJECTID DEPCODE COUNTY COUNTYNAME                DATESTAMP  ShapeSTAre ShapeSTLen
1         1      21    041  GILCHRIST 2000-05-16T00:00:00.000Z  9908353355   487300.0
2         2      54    107     PUTNAM 2000-05-16T00:00:00.000Z 23058687376   762967.7
3         3      62    123     TAYLOR 2000-05-16T00:00:00.000Z 28917470565   877252.7
4         4      46    091   OKALOOSA 2000-05-16T00:00:00.000Z 25621585036  1087058.0
5         5       7    013    CALHOUN 2000-05-16T00:00:00.000Z 16048085397   631344.0
6         6      50    099 PALM BEACH 2000-05-16T00:00:00.000Z 61491591969  1522726.9
7         7      56    111  ST. LUCIE 2000-05-16T00:00:00.000Z 16149484821   829518.7
8         8      51    101      PASCO 2000-05-16T00:00:00.000Z 21260807286   739390.6
9         9      20    039    GADSDEN 2000-05-16T00:00:00.000Z 14744963472   673612.4
10       10      37    073       LEON 2000-05-16T00:00:00.000Z 19584016138   807989.7
                         geometry
1  MULTIPOLYGON (((-82.65814 2...
2  MULTIPOLYGON (((-81.58084 2...
3  MULTIPOLYGON (((-83.73037 3...
4  MULTIPOLYGON (((-86.39159 3...
5  MULTIPOLYGON (((-84.93266 3...
6  MULTIPOLYGON (((-80.1119 26...
7  MULTIPOLYGON (((-80.23992 2...
8  MULTIPOLYGON (((-82.05534 2...
9  MULTIPOLYGON (((-84.28258 3...
10 MULTIPOLYGON (((-84.00743 3...

The object has the familar tibble-style output that we are used to (e.g. it only prints the first 10 rows of the data). However, it also has some additional information in the header, like a description of the geometry type (MULTIPOLYGON). Notice the geometry column right at the end of the data frame. This geometry column is how the sf package achieves much of its magic: It stores the geometries of each row element in its own list column. For example, we could print out the coordinates needed to plot the first element in our data frame, by typing fla_counties$geometry[[1]]. In contrast, compare with how complicated the structure of a traditional spatial object is by running, say, str(as(fla, "Spatial")).]

# this may print long output
# str(as(fla, "Spatial"))

Plotting and projection with ggplot2

Plotting sf objects is incredibly easy thanks to the package’s integration with ggplot2. Plotting sf maps using the base R plot function is also easy (and often faster). The key geom to remember is geom_sf(). For example:

# library(tidyverse) ## Already loaded
ggplot(data = fla_counties) +
  geom_sf(aes(fill = ShapeSTAre), alpha=0.8, col="white") +
  scale_fill_viridis_c(name = "Area") +
  ggtitle("Counties of Florida")

Largest counties in Florida:

Rank	Land Area	County
1	1,998.32 sq mi	Collier, FL
2	1,969.76 sq mi	Palm Beach, FL
3	1,897.71 sq mi	Miami-Dade, FL
4	1,797.84 sq mi	Polk, FL

fla_counties %>% 
  top_n(ShapeSTAre, n = 4)

Simple feature collection with 4 features and 7 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -82.10574 ymin: 25.14405 xmax: -80.03095 ymax: 28.36218
CRS:            4326
  OBJECTID DEPCODE COUNTY COUNTYNAME                DATESTAMP  ShapeSTAre ShapeSTLen
1        6      50    099 PALM BEACH 2000-05-16T00:00:00.000Z 61491591969    1522727
2       43      53    105       POLK 2000-05-16T00:00:00.000Z 56039729500    1273588
3       47      11    021    COLLIER 2000-05-16T00:00:00.000Z 55941206046    1903355
4       48      13    025       DADE 2000-05-16T00:00:00.000Z 55073375764    1774363
                        geometry
1 MULTIPOLYGON (((-80.1119 26...
2 MULTIPOLYGON (((-81.65737 2...
3 MULTIPOLYGON (((-81.27162 2...
4 MULTIPOLYGON (((-80.87314 2...

Florida Lakes

Obtained from https://geodata.dep.state.fl.us/datasets/97b765ff2b70400d8bcab23fbe2a5e88_0

loc_lakes <- "./Florida_Lakes/Florida_Lakes.shp"
## Read the shapefile into R
fla_lakes <- st_read(loc_lakes, quiet = TRUE)

Plot lakes in Polk county

dim(fla_lakes)

[1] 4243    7

colnames(fla_lakes)

[1] "PERIMETER" "NAME"      "COUNTY"    "OBJECTID"  "SHAPEAREA" "SHAPELEN"  "geometry" 

fla_lakes %>% 
  filter(COUNTY == "POLK") %>% 
  ggplot() +
  geom_sf(aes(fill = SHAPEAREA), alpha=0.8, col="white") +
  #scale_fill_viridis_c(name = "Area") +
  ggtitle("Lakes in Polk County")

Biggest lakes in Polk County

fla_lakes %>% 
  filter(COUNTY == "POLK") %>% 
  top_n(SHAPEAREA, n = 10) 

Simple feature collection with 10 features and 6 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.951 ymin: 27.66273 xmax: -81.36433 ymax: 28.10574
CRS:            4326
   PERIMETER             NAME COUNTY OBJECTID SHAPEAREA SHAPELEN                       geometry
1   16855.88    Lake Hamilton   POLK       67   8668138 16855.88 MULTIPOLYGON (((-81.65413 2...
2   16652.12       Reedy Lake   POLK       76  14191962 16652.12 MULTIPOLYGON (((-81.52262 2...
3   21815.34      Lake Parker   POLK     1608   8533253 21815.34 MULTIPOLYGON (((-81.93106 2...
4   21007.95      Lake Marion   POLK     1936  12226186 21007.95 MULTIPOLYGON (((-81.53609 2...
5   31582.90      Lake Pierce   POLK     1937  15604508 31582.90 MULTIPOLYGON (((-81.49978 2...
6   34664.94     Crooked Lake   POLK     2690  16885793 34664.94 MULTIPOLYGON (((-81.57241 2...
7   21529.93 Lake Weohyakapka   POLK     2700  30549713 21529.93 MULTIPOLYGON (((-81.4103 27...
8   18557.39     Lake Rosalie   POLK     2701  18358611 18557.39 MULTIPOLYGON (((-81.39447 2...
9   19950.41     Lake Hancock   POLK     2703  18549186 19950.40 MULTIPOLYGON (((-81.83617 2...
10  20955.50    Lake Arbuckle   POLK     3907  15412988 20955.50 MULTIPOLYGON (((-81.41045 2...

To reproject an sf object to a different CRS, we can use sf::st_transform().

fla_counties %>%
  st_transform(crs = "+proj=moll +ellps=WGS84") %>% ## Reprojecting to a Mollweide CRS
  head(2) ## Saving vertical space

Simple feature collection with 2 features and 7 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -7625474 ymin: 3564689 xmax: -7485526 ymax: 3636317
CRS:            +proj=moll +ellps=WGS84
  OBJECTID DEPCODE COUNTY COUNTYNAME                DATESTAMP  ShapeSTAre ShapeSTLen
1        1      21    041  GILCHRIST 2000-05-16T00:00:00.000Z  9908353355   487300.0
2        2      54    107     PUTNAM 2000-05-16T00:00:00.000Z 23058687376   762967.7
                        geometry
1 MULTIPOLYGON (((-7586204 36...
2 MULTIPOLYGON (((-7486900 36...

Or, we can specify a common projection directly in the ggplot() call using coord_sf(). This is often the most convenient approach when you are combining multiple sf data frames in the same plot.

ggplot(data = fla_counties) +
  geom_sf(aes(fill = ShapeSTAre), alpha=0.8, col="white") +
  scale_fill_viridis_c(name = "Area") +
  coord_sf(crs = "+proj=moll +ellps=WGS84") +
  labs(
    title = "Counties of Florida",
    subtitle = "Mollweide projection"
    ) 

Data wrangling with dplyr and tidyr

The tidyverse approach to data wrangling carries over very smoothly to sf objects. For example, the standard dplyr verbs like filter(), mutate() and select() all work. You can also perform group_by() and summarise() operations as usual. Furthermore, the dplyr family of join functions also work, which can be especially handy when combining different datasets by (say) Federal Information Processing Standard (FIPS) code or some other attribute. However, this presumes that only one of the objects has a specialized geometry column. In other words, it works when you are joining an sf object with a normal data frame. In cases where you want to join two sf objects based on their geometries, there is a specialized st_join() function.

Specialized geometric operations

Alongside all the tidyverse functionality, the sf package comes with a full suite of geometrical operations. You should take a look at the third sf vignette or the Geocomputation with R book to get a complete overview.

So-called unary operations are applied to a single object. For instance, you can “melt” sub-elements of an sf object (e.g. counties) into larger elements (e.g. states) using sf::st_union():

fla_counties %>% 
  st_union() %>% 
  ggplot() +
  geom_sf(fill=NA, col="black") +
  labs(title = "Outline of Florida") 

Or, you can get the st_area(), st_centroid(), st_boundary(), st_buffer(), etc. of an object using the appropriate command. For example:

fla_counties %>% st_area() %>% head(5) ## Only show the area of the first five counties to save space.

Units: [m^2]
[1]  920490582 2142439339 2685252824 2369242717 1487627072

And:

fla_centroid <- st_centroid(fla_counties)

ggplot(fla_counties) +
  geom_sf(fill = "black", alpha = 0.8, col = "white") +
  geom_sf(data = fla_centroid, col = "red") + ## Notice how easy it is to combine different sf objects
  labs(
    title = "Counties of Florida",
    subtitle = "Centroids in red"
    )

Another set of so-called binary operations can be applied to multiple objects. So, we can get things like the distance between two spatial objects using sf::st_distance(). In the below example, we get the distance from Polk county to Palm Beach county:

fla_counties %>% 
  filter(COUNTYNAME %in% c("POLK", "PALM BEACH"))

Simple feature collection with 2 features and 7 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -82.10574 ymin: 26.32093 xmax: -80.03095 ymax: 28.36218
CRS:            4326
  OBJECTID DEPCODE COUNTY COUNTYNAME                DATESTAMP  ShapeSTAre ShapeSTLen
1        6      50    099 PALM BEACH 2000-05-16T00:00:00.000Z 61491591969    1522727
2       43      53    105       POLK 2000-05-16T00:00:00.000Z 56039729500    1273588
                        geometry
1 MULTIPOLYGON (((-80.1119 26...
2 MULTIPOLYGON (((-81.65737 2...

polk_palmbeach <- fla_counties %>% 
                    filter(COUNTYNAME %in% c("POLK", "PALM BEACH"))
palmbeach <- fla_counties %>% 
                filter(COUNTYNAME %in% c("PALM BEACH"))

## Use "by_element = T" to give a vector instead of the default pairwise matrix
ab_dist <- st_distance(polk_palmbeach, palmbeach, by_element = T)
     

## We can use the `units` package (already installed as sf dependency) to convert to kilometres 
ab_dist <- ab_dist %>% units::set_units(km) %>% round()

Create plot:

ggplot(fla_counties) +
  geom_sf(fill = "black", alpha = 0.8, col = "white") +
  geom_sf(data = fla_counties %>% filter(COUNTYNAME %in% c("POLK", "PALM BEACH")), 
          aes(fill = COUNTYNAME), col = "white") +  
  labs(
    title = "Calculating distances",
    subtitle = paste0("The distance between Polk and Palm Beach is ", ab_dist[2], " km")
    ) +
  theme(legend.title = element_blank())

Getting Map Data

We can easily import external shapefiles, KML files, etc., into R using the generic sf::st_read() function.

However, R also provides access to a large number of base maps — e.g. countries of the world, US states and counties, etc. — through the maps, (higher resolution) mapdata and spData packages, as well as a whole ecosystem of more specialized GIS libraries.

The World

# library(maps) ## Already loaded
world1 <- st_as_sf(map("world", plot = FALSE, fill = TRUE))
world1 %>%
  ggplot() + 
  geom_sf(fill = "grey80", col = "grey40", lwd = 0.3) +
  labs(title = "The world", subtitle = "Mercator projection")

All of the usual sf functions and transformations can then be applied. For example, we can reproject the above world map onto the Lambert Azimuthal Equal Area projection (and further orientate it at the South Pole) as follows.

world2 <- 
  st_transform(
  world1,
  "+proj=laea +y_0=0 +lon_0=155 +lat_0=-90 +ellps=WGS84 +no_defs"
  )

world2 %>%
  ggplot() + 
  geom_sf(fill = "grey80", col = "grey40", lwd = 0.3) +
  labs(title = "The world", subtitle = "Lambert Azimuthal Equal Area projection")

A digression on the Winkel tripel projection

As we have already seen, most map projections work great “out of the box” with sf. One notable exception is the Winkel tripel projection. This is the preferred global map projection of National Geographic and requires a bit more work to get it to play nicely with sf and ggplot2 (as detailed in this thread). Here’s a quick example of how to do it:

# library(lwgeom) ## Already loaded

world3 <- lwgeom::st_transform_proj(world1, crs = "+proj=wintri +datum=WGS84 +no_defs +over")
## Don't necessarily need a graticule, but if you do then define it manually:
gr <- 
  st_graticule(lat = c(-89.9,seq(-80,80,20),89.9)) %>%
  lwgeom::st_transform_proj(crs = "+proj=wintri +datum=WGS84 +no_defs +over")

world3 %>%
  ggplot() + 
  geom_sf(data = gr, color = "#cccccc", size = 0.15) + ## Manual graticule
  geom_sf(fill = "grey80", col = "grey40", lwd = 0.3) +
  coord_sf(datum = NA) +
  labs(title = "The world", subtitle = "Winkel tripel projection")

A single country (e.g. Colombia)

The maps and mapdata packages have detailed county- and province-level data for several individual nations. We can still use it to extract a specific country’s border using some intuitive syntax. For example, we could plot a base map of Colombia as follows.

colombia <- st_as_sf(map("world", "colombia", plot = FALSE, fill = TRUE))

colombia %>%
  ggplot() + 
  geom_sf(fill="blue", col=NA)

Aside: You can detach the maps package once you are finished using it, since it avoids potential namespace conflicts with purrr::map.

detach(package:maps) ## To avoid potential purrr::map() conflicts

US Census data with tigris and tidycensus

Working with Census data may require some extra work. You need to register on the website, then download data from various years or geographies separately, merge these individual files, etc. Thankfully, this has recently become much easier thanks to the Census API and — for R at least — the tigris and tidycensus packages from Kyle Walker. Check his tutorial on his website.

NOTE: Before continuing with this section, you will first need to request an API key from the Census. Once that’s done, you can set it using the tidycensus::census_api_key() function. I recommend using the “install = T” option to save your key for future use.

We start by loading the packages and setting some options for optimized use with the sf package.

# library(tidycensus) ## Already loaded
# library(tigris) ## Already loaded
options(tigris_class = "sf")
options(tigris_use_cache = TRUE)
# census_api_key("YOUR KEY HERE", install = T)

Let’s say that our goal is to provide a snapshot of Census rental estimates across different cities in Florida. We start by downloading tract-level rental data for Florida using the tidycensus::get_acs() function. Note that you will need to look up the correct ID variable (in this case: “DP04_0134”).

rent <- 
  tidycensus::get_acs(
    geography = "tract", variables = "DP04_0134",
    state = c("FL"), geometry = TRUE
    )

Getting data from the 2014-2018 5-year ACS
Using the ACS Data Profile
Using FIPS code '12' for state 'FL'

rent

Simple feature collection with 4245 features and 5 fields (with 33 geometries empty)
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -87.63494 ymin: 24.5231 xmax: -80.03136 ymax: 31.00089
CRS:            4269
First 10 features:
         GEOID                                       NAME  variable estimate moe
1  12001000200    Census Tract 2, Alachua County, Florida DP04_0134      853  74
2  12001000301 Census Tract 3.01, Alachua County, Florida DP04_0134      840 105
3  12001000302 Census Tract 3.02, Alachua County, Florida DP04_0134      933  92
4  12001000400    Census Tract 4, Alachua County, Florida DP04_0134     1019 432
5  12001000500    Census Tract 5, Alachua County, Florida DP04_0134      861  90
6  12001000600    Census Tract 6, Alachua County, Florida DP04_0134      865 169
7  12001000700    Census Tract 7, Alachua County, Florida DP04_0134      835  41
8  12001000806 Census Tract 8.06, Alachua County, Florida DP04_0134      786  54
9  12001000808 Census Tract 8.08, Alachua County, Florida DP04_0134      906  43
10 12001000809 Census Tract 8.09, Alachua County, Florida DP04_0134      893  48
                         geometry
1  MULTIPOLYGON (((-82.33935 2...
2  MULTIPOLYGON (((-82.33922 2...
3  MULTIPOLYGON (((-82.33916 2...
4  MULTIPOLYGON (((-82.32413 2...
5  MULTIPOLYGON (((-82.33082 2...
6  MULTIPOLYGON (((-82.31147 2...
7  MULTIPOLYGON (((-82.32629 2...
8  MULTIPOLYGON (((-82.33938 2...
9  MULTIPOLYGON (((-82.34001 2...
10 MULTIPOLYGON (((-82.35586 2...

This returns an sf object, which we can plot directly.

rent %>%
  ggplot() + 
  geom_sf(aes(fill = estimate, color = estimate)) + 
  scale_fill_viridis_c(name = "Rent ($)", labels = scales::comma) + 
  scale_color_viridis_c(name = "Rent ($)", labels = scales::comma) +
  labs(
    title = "Rental rates across Florida", 
    caption = "Data: US Census Bureau"
    ) 

The above map is very detailed. Perhaps we are not interested in the entire set of tract level data, but would rather get a sense of rents within some well-defined metropolitan areas? The tigris package has you covered here. For example, let’s say we want to compare average rents across three metros: Orlando, Tampa and Miami.

fl_metros <- 
  tigris::core_based_statistical_areas(cb = TRUE) %>%
  # filter(GEOID %in% c("53000", "71000", "45000")) %>% ## Could use GEOIDs directly if you know them 
  filter(grepl("Orlando|Tampa|Miami", NAME)) %>%
  select(metro_name = NAME)

Now we do a spatial join on our two data sets using the sf::st_join() function.

fl_rent <- 
  st_join(
    rent, 
    fl_metros, 
    join = st_within, left = FALSE
    ) 
fl_rent

Simple feature collection with 2349 features and 6 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -82.85243 ymin: 25.13807 xmax: -80.03136 ymax: 29.27677
CRS:            4269
First 10 features:
          GEOID                                         NAME  variable estimate moe
223 12011010102 Census Tract 101.02, Broward County, Florida DP04_0134     1728 393
224 12011010103 Census Tract 101.03, Broward County, Florida DP04_0134     1138 138
225 12011010104 Census Tract 101.04, Broward County, Florida DP04_0134     1341 180
226 12011010200    Census Tract 102, Broward County, Florida DP04_0134     1262  79
227 12011010304 Census Tract 103.04, Broward County, Florida DP04_0134     1074 199
228 12011010305 Census Tract 103.05, Broward County, Florida DP04_0134     1146  45
229 12011010306 Census Tract 103.06, Broward County, Florida DP04_0134     1042  98
230 12011010307 Census Tract 103.07, Broward County, Florida DP04_0134     1142  59
231 12011010308 Census Tract 103.08, Broward County, Florida DP04_0134     1940 119
232 12011010401 Census Tract 104.01, Broward County, Florida DP04_0134     1732 102
                                   metro_name                       geometry
223 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.09418 2...
224 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.09123 2...
225 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.08109 2...
226 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.10748 2...
227 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.11944 2...
228 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.12199 2...
229 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.10789 2...
230 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.12833 2...
231 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.15345 2...
232 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.17015 2...

One useful way to summarize this data and compare across metros is with a histogram. Note that “regular” ggplot2 geoms and functions play perfectly nicely with sf objects (i.e. we are not limited to geom_sf()).

fl_rent %>%
  ggplot(aes(x = estimate)) + 
  geom_histogram() + 
  facet_wrap(~metro_name, labeller = label_wrap_gen(width=10)) + 
  theme_bw()